Simulational nanoengineering: Molecular dynamics implementation of an atomistic 

Stirling engine 
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A nanoscale-sized Stirling engine with an atomistic working fluid has been modeled using molec- 
ular dynamics simulation. The design includes heat exchangers based on thermostats, pistons at- 
tached to a flywheel under load, and a regenerator. Key aspects of the behavior, including the 
time-dependent flows, are described. The model is shown to be capable of stable operation while 
producing net work at a moderate level of efficiency. 

PACS numbers: 05.70.Ln, 47.61.-k, 02.70.Ns 



The Stirling engine, an external combustion engine in- 
vented almost two centuries ago, and an early competitor 
of the steam engine, continues to attract interest owing 
to its potential effectiveness as a power source and its 
success in specialized applications [H, S Q . The under- 
lying thermodynamic cycle consists of four stages: low- 
temperature isothermal compression of the working fluid, 
constant-volume displacement of the fluid between the 
cold and hot spaces of the engine with (optionally) the 
fluid gaining heat while passing through a heat-storing 
regenerator, high-temperature isothermal expansion, and 
constant- volume displacement between the hot and cold 
spaces with heat being returned to the regenerator. The 
net work is the excess energy produced by the expand- 
ing hot fluid over that needed to recompress the cold 
fluid. Unlike internal combustion engines, there is com- 
plete flexibility as to fuel type, including even solar en- 
ergy; the fact that the complexities of combustion 0] are 
avoided makes the Stirling engine an attractive candidate 
for simulation. 

This paper describes the molecular dynamics (MD) 
simulation of a simplified Stirling engine using a discrete- 
particle working fluid. MD [f| provides the capabil- 
ity for direct atomistic modeling of nanomachinery, to- 
gether with the accompanying complex thermodynamic 
and fluid-dynamic processes. While past MD studies 
have included instances of heat production as a byprod- 
uct of work, e. g., i n fracture |8|, Taylor-Couette flow [9( 
and fluid jets the present simulations extend the 

scope of MD to systems that harness thermal energy for 
producing useful workfl2l|. 

The model engine incorporates a number of features: 
pistons driven by collisions with the fluid atoms, heat 
input and output controlled by thermostats that main- 
tain isothermal conditions in the hot and cold spaces, 
a rotating flywheel connected to the pistons and sub- 
ject to an applied load, and an atomistic regenerator in- 
tended to reduce heat wastage. The working fluid ex- 
periences substantial temperature and pressure variation 
over the cycle, but its thermodynamic and heat/mass 
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FIG. 1: (Color online) Model Stirling engine, showing hot (left 
side) and cold spaces, regenerator (top center), pistons, coun- 
terclockwise rotating flywheel, linkages, and working fluid. 



transport properties are intrinsic to the atomistic model. 
This represents an advantage of MD over the alternative 
continuum-based analysis where such details must 
be provided separately; if these and subsequent more ex- 
tensive studies can be quantitatively validated, the role 
of MD in nanoengineering is assured. 

There are different Stirling engine designs, each in- 
volving engineering considerations that are not addressed 
here. It is the two-piston alpha version, shown in Fig. [T] 
that is modeled (the beta version, with a single piston and 
a mechanical displacer, is one of the alternatives). Even 
though the simulations consider a system that is an ideal- 
ization of a real Stirling engine, its relative complexity - 
by MD standards - calls for several implementation deci- 
sions, with the goal of achieving a simple but reasonably 
complete operational model. The two-dimensional case 
is considered for convenience. 

The atoms of the working fluid are soft-disk parti- 
cles, interacting through a short-range potential U(r) — 
4e[(cr/r) 12 - (cr/r) 6 ] + e with r < r c = 2 1 / 6 er (condensa- 
tion is avoided since there is no attraction); similar disks 
line the fixed engine walls and the moving piston faces. 
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The linkages connecting the pistons to the flywheel en- 
sure sinusoidal piston motion, but are idealized since the 
transverse parts are of variable length; the phase angle 
between the pistons is fixed at 90°. The fluid atoms and 
movable engine components obey Newtonian dynamics; 
the force evaluation and solution of the equations of mo- 
tion employ standard MD techniques [7(. 

In reduced MD units, the model engine has stroke 
2Rf = 112.4, where Rf is the flywheel radius and bore 
67.7 (in MD units where, typically, a = 1 corresponds to 
a length unit of 3.4A, e = 1 to a time unit of 2.16 psec if 
fluid atoms have mass m = 1, and ks — 1 defines temper- 
ature) ; the actual value of Rf amounts to a mere 19lA. 
The working fluid contains N a — 1368 atoms; this 'large' 
system is the main focus of the analysis, but a half-size 
'small' system (shown in Fig. []} with N a — 407 is also 
considered. The minimal size enables runs of adequate 
length without excessive computation, while demonstrat- 
ing the extreme scales at which behavior remains reason- 
able. 

Heat exchangers responsible for thermal input and out- 
put arc modeled using thermostats. These act sepa- 
rately over the hot and cold spaces swept by the pis- 
tons, are applied every time step, and entail rescaling 
the velocities of all the fluid atoms in each space, af- 
ter allowing for the current average flow. Use of ther- 
mostats, rather than thermalizing walls, prevents tem- 
perature fluctuations and enhances heat transfer (since 
unrestricted transfer rates from and to exogenous heat 
reservoirs are implied). 

The flywheel phase angle 8 is defined so that the hot 
piston is at its top position when 8 = 0°. Average fluid 
density N a /A(9), where A(8) is the area occupied by the 
fluid, varies between 0.256 and 0.085, corresponding to a 
3.03:1 compression ratio. The inertial contribution of the 
drive mechanism is assigned to the flywheel, with pistons 
and linkages assumed massless. The flywheel moment 
of inertia // = (Tr/2)pfR*; for density pf — 10 4 the 
rotation speed u> is low enough that piston speeds (< 
RfUi) are typically an order of magnitude less than the 
mean thermal velocity of the fluid atoms vt (the resulting 
slow fluid flow justifies the thermostat implementation). 
The flywheel kinetic energy also exceeds that of the fluid 
by a similar factor, ensuring minimal fluctuations in to. 

The role of the regenerator was recognized by its inclu- 
sion in the early Stirling design. Although not essential 
for operation, its function is to improve heat utilization 
by extracting heat from the hot fluid as it flows to the 
cold space and returning it when flow is reversed; typi- 
cally implemented as a metallic mesh, it also impedes the 
flow. The model regenerator is assembled from N r = 32 
atoms with m — 50, tethered to a uniform grid by a 
potential 50U(r c — r), where r is the distance from the 
relevant grid site. Transiting fluid atoms collide with the 
(unthermostatted) tethered atoms and heat is transferred 
to or from the regenerator, depending on the relative 
temperatures; the heavier atoms undergo limited, rela- 
tively slow motion, but the contribution to performance 
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FIG. 2: Average rotation speed u> (mrad/time), for runs at 
different temperatures T^; N c is the number of cycles. 



0.14 




0.04 

123456789 10 

FIG. 3: Dependence of relative efficiency e/eo on Th for large 
and small systems. 



in terms of heat storage is minimal since, to lessen flow 
impedance, N r <C N a . The area occupied by the regener- 
ator is connected to the hot and cold spaces by triangular 
ducts. 

Work is performed against an applied torque load 
jRfUJ, with 7 = 100. There is no energy loss due to 
the rough boundary walls since all collisions are elastic 
(energy is conserved when thermostats and load are re- 
moved). Temperature in the hot space Th is a parame- 
ter, and the cold space T c is fixed at unity. The system 
is started at maximum compression with the pistons at 
equal height (8 — 45°) and rotation speed ujq. Run length 
is 2 x 10 9 time steps (of size 0.005), so while all runs cover 
the same time period, the number of cycles differs. 

The principal outcome of the present simulations is a 
demonstration of the operational capability of the en- 
gine. Fig. [5] shows the cycle-averaged rotation speed Q 
as a function of the number of cycles N c , for runs at 
different Th (in real units, typical speeds are ~ 5 x 10 7 
rot /s) . lu converges to a steady value following a transient 
stage lasting < 200 cycles, with no long-term drift and 
only small size-dependent fluctuations. In addition to Q 
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FIG. 4: Deviation from average rotation speed Slu as a func- 
tion of phase angle 6. 



FIG. 6: Regenerator temperature T r (ff). 




FIG. 5: Heat input and output, Qi{6) and —Q o {0) (the latter 
negative for clarity), over the cycle. 
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FIG. 7: Pressure P(9) as a function of fluid area A(6)/N a , 
with and without regenerator (Tj» = 4) . 



increasing with T/j, it also depends on the other parame- 
ters: for example, a higher load (7), but not too high to 
prevent rotation, reduces uj; the same is true if regenera- 
tor collisional flow impedance is increased by raising N r ; 
uj increases with fluid density (by varying N a ) over a rea- 
sonable range. The choice of ojq (even uiq = 0) does not 
affect the outcome, provided the flywheel is able to com- 
plete the initial rotations, otherwise there are damped 
oscillations about the state of maximum expansion. The 
overall behavior of the small system is similar, though 
with enhanced fluctuations. 

Engine performance is expressed in terms of thermody- 
namic efficiency, the ratio of net work W to heat input Qi, 
the latter represents the true operational cost, whereas 
heat output to the environment Q is waste. Qi and Q Q 
are the measured kinetic energy changes produced by the 
thermostats, and the source of W is the resistive torque 
on the flywheel; since there are no other loss mechanisms, 
W = Qi — Q , where the averaging (over the entire run, 
excluding the initial transient stage) allows for fluctua- 
tions. Efficiency e is expressed relative to the thermo- 
dynamic (Carnot) limiting value cq — 1 — T c /Th- Fig. [3] 



shows the relative efficiency e/eo as a function of 7},, for 
both large and small systems (size dependence is signif- 
icant at such scales), with values reaching a respectable 
0.13 (13%). The presence of the regenerator degrades 
performance, and when removed e/eo is increased to 0.24 
at Th =4; higher 7 also raises e since more work is done. 

The dependence of rotation speed (run-averaged) on 9 
is shown in Fig. [4j since the variation is weak (greater, 
even twofold, variation occurs in the small system), the 
deviation 5lo{9) = lo{9) — uj is shown. Speed is nonsi- 
nusoidal, slowest at 9 ~ 45° at maximum overall com- 
pression, and fastest at 9 ~ 200° approaching maxi- 
mum expansion; during the power portion of the cycle, 
duj / d9 > when work performed by the engine exceeds 
load. The range of Sui is reduced at higher T^; likewise 
if If is increased. 

Further measurements of 9 dependence provide a more 
detailed picture of the behavior. Fig. [5] shows heat input 
and output through the thermostats, Qi{9) and Q o (0)', 
due to significant fluctuations, averaging over multiple 
cycles is again necessary. The details depend on T^, with 
maximal Qi(9) at 6 ~ 90°, close to where |Q O (0)| is at 
a minimum, and a value near zero at 9 ~ 270°; both 





FIG. 8: (Color online) Fluid flow during the cycle (from upper left, 9 = 0°, 45°,. . . ); to resolve the wide range of flow speeds, 
arrow size oc (speed) 1//2 (slower flows appear faster). Flow features are discussed in the text. 



also exhibit small irregularities in the vicinity of 9 ~ 0° 
where the fluid is compressed and the hot space practi- 
cally empty. 

The effective regenerator temperature T r (9) is derived 
from the average kinetic energy of the tethered atoms; 
the variation over the cycle is shown in Fig. [S] For 
Th > 2 there is a minimum at 9 ~ 75° as the cold pis- 
ton approaches the top position, and a broader maximum 
around 9 ~ 270°; both limiting values increase with Th- 
Small irregularities appear near 0°, as before, that are 
consequences of the model design and more pronounced 
in the small system. The temperature variation shows 
that the regenerator responds correctly, although the ca- 
pacity for heat storage is very limited. Heat throughput 
is increased (together with e) when the regenerator is 
removed to lower flow impedance. 

Fluid pressure is evaluated from the virial; Fig. 
shows pressure variation over the cycle in the hot and 
cold spaces, Ph{9) and P c (9), as a function of fluid area 
A(9)/N a , for runs at Th = 4 with and without the re- 
generator; although similar to an indicator diagram Q, 
the area enclosed does not correspond to W away from 
equilibrium. The general absence of pressure equalization 



between hot and cold spaces is primarily due to regener- 
ator flow impedance, with p^ax/pmin ~ 2 P™ ax /P™ in - 
without the regenerator Ph(9) ~ P c {9). Pressure depen- 
dence on 9 (not shown) is not sinusoidal, with a range 
of variation that increases with Th] maxima are at 9 be- 
tween 50° and 80°, in advance of the u>(9) inflection point 
at ~ 120° where power output begins to decline, with 

prnax prece di n g P™ ax by ~ 10°. 

Finally, at the most detailed level of examination, the 
cyclically changing fluid flow is evaluated from spatially 
coarse-grained velocity fields accumulated for 9 segments 
of range 15° and averaged over 100 successive cycles to 
reduce noise (some small variation persists even after av- 
eraging). The frames of Fig. [5] show every third segment 
(at Th = 4). Since bulk flow rates are smaller than vt by 
at least an order of magnitude, certain flow features only 
become apparent after averaging: because of thermal ex- 
pansion, at 135° (frame 4), fluid can be seen leaving the 
hot space and flowing back through the regenerator, even 
though the hot piston is still descending, with the reverse 
occurring at 315° (frame 8); there is also a vortex that 
forms near the top of the inner wall of the cold space 
while the cold piston is descending (frames 5-7). 
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